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A Survey of Trefftz Methods for the Helmholtz 
Equation 
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Abstract Trefftz methods are finite element-type schemes whose test and trial func¬ 
tions are (locally) solutions of the targeted differential equation. They are particu¬ 
larly popular for time-harmonic wave problems, as their trial spaces contain oscillat¬ 
ing basis functions and may achieve better approximation properties than classical 
piecewise-polynomial spaces. 

We review the construction and properties of several Trefftz variational formula¬ 
tions developed for the Helmholtz equation, including least squares, discontinuous 
Galerkin, ultra weak variational formulation, variational theory of complex rays and 
wave based methods. The most common discrete Trefftz spaces used for this equa¬ 
tion employ generalised harmonic polynomials (circular and spherical waves), plane 
and evanescent waves, fundamental solutions and multipoles as basis functions; we 
describe theoretical and computational aspects of these spaces, focusing in particu¬ 
lar on their approximation properties. 

One of the most promising, but not yet well developed, features of Trefftz meth¬ 
ods is the use of adaptivity in the choice of the propagation directions for the basis 
functions. The main difficulties encountered in the implementation are the assem¬ 
bly and the ill-conditioning of linear systems, we briefly survey some strategies that 
have been proposed to cope with these problems. 
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1 Introduction 

Given a linear PDE, a Trefftz method is a volume-oriented discretisation scheme, for 
which all trial and test functions, when restricted to any element of a given mesh, are 
solutions of the PDE under consideration. The name comes from the work Oa of 
E. Trefftz, dating back to 1926, where this idea was applied to the Laplace equation. 
Since then, several versions of Trefftz methods have been proposed and applied to a 
range of PDEs by different groups of mathematicians, engineers and computational 
scientists, often unaware of each other. Typical PDEs addressed are linear, with 
piecewise-constant coefficients and homogeneous, i.e. with vanishing volume source 
term. 

Trefftz methods are related to both hnite element (EEM) and boundary element 
methods (BEM). With the former they have in common that they provide a dis¬ 
cretisation in the volume. With the latter they share some characteristics such as the 
need of integration on lower-dimensional manifolds only. Compared to conventional 
EEMs, Trefftz methods have attracted attention mainly for two reasons: (i) they may 
need much fewer degrees of freedom than standard schemes to achieve the same ac¬ 
curacy, and (ii) they incorporate some properties of the problem’s solution (such as 
oscillatory character, wavelength, maximum principle, boundary layers) in the trial 
spaces, and thus also in the discrete solution. In addition, compared to BEMs, an 
advantage of Trefftz schemes is that they do not require the evaluation of singular 
integrals. 

Comparing with finite and boundary elements, in 1997 Zienkiewicz II121I stated: 
“... it seems without doubt that in the future Trefftz type elements will frequently 
be encountered in general finite element codes.... It is the author’s belief that the 
simple Trefftz approach will in the future displace much of the boundary type anal¬ 
ysis with singular kernels.” While this prediction has not yet come true, in the last 
years more and more work has been devoted to the formulation, the analysis and the 
validation of these methods and substantial progress has been accomplished. 

In this chapter we survey Trefftz hnite element methods for the homogeneous 
Helmholtz equation (—Au — k^u = 0), which models acoustic wave propagation in 
time-harmonic regime. Eor medium and high frequencies, i.e. for values of kL in a 
range of 10^ to 10^, where k > 0 is the wavenumber, and L a characteristic length of 
the region of interest, the numerical solution of the Helmholtz equation in 2D and 
3D is particularly challenging. A main reason is that Helmholtz solutions oscillate 
with a wavelength proportional to the inverse of k. Hence, piecewise polynomials do 
not provide efficient approximation. Trefftz schemes are thus particularly relevant 
as they can improve on the point where (polynomial) EEMs fail: the approximation 
properties of the basis functions. Moreover, some Trefftz methods can remedy other 
shortcomings that often haunt discretisations of time-harmonic problems, such as 
the lack of coercivity and the presence of minimal resolution conditions to guarantee 
unique solvability. Theorem |2] in this chapter is an example. Earlier overviews of 
Trefftz schemes for the Helmholtz equation, together with numerous references, 
can be found in ll98l , lIFSl Ch. 1] and ll7^ Ch. 3]. Surveys of Trefftz schemes for 
other equations are in 067ll75ll99lll21l . 
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For most of the Trefftz spaces used, continuity across interfaces separating mesh 
elements cannot be enforced strongly, as Trefftz functions are not as “flexible” as 
piecewise polynomials. As a consequence, the standard Helmholtz variational for¬ 
mulation posed in subspaces of the Sobolev space //* is not applicable and discreti¬ 
sations must be used that can accommodate discontinuous trial functions. A wide 
array of different variational formulations has been proposed and in ij2]we attempt 
a classification and a comparison of the best known. We identify three main classes 
of formulations; (i) least squares (LS, where squares of suitable norms of 

residuals are minimised; (ii) discontinuous Galerkin (DG, 112.21 1. whose formula¬ 
tions arise from local integration by parts and which may or may not use Lagrange 
multipliers on mesh interfaces; (Hi) weighted residual ( 112.3b . which are defined by 
testing residuals against suitable traces of test functions. The methods discussed 
include: the Trefftz-discontinuous Galerkin (TDG), the ultra weak variational for¬ 
mulation (UWVF), the discontinuous enrichment method (DEM), the variational 
theory of complex rays (VTCR) and the wave based method (WBM). Moreover, in 
the spirit of the symposium that led up to the present volume, to “build bridges” with 
a wider portion of the literature and of the computational PDE community, in 112.41 
we describe some older Trefftz schemes defined on a single element and in 112.51 
we consider some methods that are not Trefftz but use oscillating basis functions 
that are “approximately Trefftz”, such as the partition of unity method (PUM). To 
easily compare them, we write all formulations for the same Robin-Dirichlet model 
boundary value problem (see 111.1b . 

In lj2]we completely gloss over the choice of basis functions and discrete spaces 
employed, whose description is postponed to l|3] This is because, apart from few ex¬ 
ceptions such as unbounded elements, any Trefftz discrete space can be employed in 
any Trefftz variational formulation. We believe that separating the discussion of the 
two main components in the definition of a Trefftz method, i.e. variational formu¬ 
lations and discrete spaces, will make the presentation clearer. The most common 
basis functions for Trefftz methods are plane waves (x i—e’^'** * for a fixed unit vector 
d) and generalised harmonic polynomials (i.e. circular/spherical waves, products of 
circular/spherical harmonics and Bessel functions), for which quite a complete ap¬ 
proximation theory exists, see 113.1113.21 Other basis functions include fundamental 
solutions, multipoles, evanescent waves and corner waves. We note that, since the 
Helmholtz operator is the sum of a second- and a zero-order term, no non-vanishing 
piecewise-polynomial Trefftz function is possible. 

In this chapter we state a few theorems, none of them is entirely new. Lemma[T] 
exemplifies the technique of ll89l to control the norm of Trefftz functions with 
mesh-dependent norms containing interface jumps. If a Trefftz method is well-posed 
in a suitable skeleton norm, this allows to control the error in the volume; we do this 
for the LS method in Theorem [T| and for the TDG method (well-posed by Theo¬ 
rem |2]l in Corollary [T] This can be combined with the approximation results for 
circular/spherical and plane waves in il3.1Lil3.21 In brief; we provide the tools to 
derive stability and orders of L^-convergence in the volume for all Trefftz methods 
that are well-posed in suitable skeleton norms. 


4 


Ralf Hiptmair, Andrea Moiola and Ilaria Perugia 


Trefftz methods suffer from two main problems; ill-conditioning due to the poor 
linear independence of the basis functions, and the need for numerical quadrature 
for oscillating integrands. On the other hand, since the PDE is solved exactly in each 
element, only low-dimensional integrals on the mesh skeleton need to be evaluated, 
leading to massively reduced computational cost for the assembly of the linear sys¬ 
tems. Moreover, if plane wave bases are used, on any polygonal/polyhedral mesh 
the integrals can be computed analytically in a cheap way. In ij4]we briefly review 
strategies developed to deal with the computation of matrix entries and to cope with 
ill-conditioning. 

Some Trefftz methods also provide an attractive framework for implementing 
non-standard adaptive policies, like directional adaptivity following dominant wave 
directions. This is made possible, because plane wave-type Trefftz functions natu¬ 
rally encode a direction of propagation. More details are given in ^4.21 

As mentioned, in this chapter we only discuss the Helmholtz equation, i.e. acous¬ 
tic problems, and constant material parameters. The discrete Trefftz spaces used 
for the Helmholtz equation with variable coefficients are briefly addressed in 113.41 
Other time-harmonic wave problems that have been tackled with Trefftz methods in¬ 
clude electromagnetism (Maxwell equations) 11811851 . linearised Euler equation and 
general hyperbolic systems (SI, linear elasticity (Navier equation) m, (fourth or¬ 
der) Kirchhoff-Love plates 02711701176)11001 . Koiter’s linear shell theory 01001 . poro- 
elasticity IEtI §5.4], coupled vibro-acoustic problems IEt). a list of applications 
and references can be found in 024) §5.1] (with a focus in vibrational mechanics) 
and in 076II85I . A related application is tackled by the method of particular solu¬ 
tions (MPS) of 1)161 [^ . which uses Helmholtz solutions to approximate Laplace 
eigenvalue problems; in this setting the wavenumber is part of the unknowns. Eor 
recent work on space-time Trefftz methods for wave propagation in time-domain 
see ll69l and references therein. 

Several comparisons of the numerical performances of different Trefftz schemes 
for simple model problems have been published, e.g. Q (PUM, DEM, gener¬ 
alised EEM), gOl (LS, UWVE), JgO] (PUM, UWVE), l|39l (DG, UWVE, LS), HB] 
(DEM, UWVE, PUM), ||59l (LS, UWVE, VTCR), where we have included the PUM 
even if strictly speaking it is not a Trefftz method. However, from these results it is 
difficult to conclude that any formulation is clearly preferable from a computational 
point of view. A general conclusion might be that, in order to achieve the best ac¬ 
curacy and conditioning, the choice of the approximation space matters more than 
that of the variational formulation. We reiterate that these two choices are mutu¬ 
ally independent; any Trefftz discrete space might be used in any Trefftz variational 
formulation. We make some further concluding remarks in Sj^ 


1.1 Model boundary value problem 

We rely on a simple model boundary value problem (BVP) for the Helmholtz equa¬ 
tion that will be used to describe and compare the different Trefftz methods. Let 









A SuiTey of Trefftz Methods for the Helmholtz Equation 


5 


Q C K", n = 2,3, be a bounded, Lipschitz, connected domain, with dO. = FoUrj}, 
where Fd and Fr are disjoint components of dO.; Fr^& while Fd might be empty. 
Denote by n the outward-pointing unit normal vector field on dQ. We consider the 
homogeneous Robin-Dirichlet BVP 

—Au —k^u — 0 

U=gD 
du 

— + ikT^U = gR 

dn 

Here go and gR are the boundary data, i is the imaginary unit, k (the wavenum¬ 
ber) and d (the impedance parameter) are positive constants. We assume that Q, go 
and gR are such that u G for some i > 0. In typical sound-soft acous¬ 

tic scattering problems, Fo represents the boundary of the scatterer, and Fr stands 
for an artificial truncation of the unbounded region where waves propagate; see 

e.g. m §2]. 

Simple generalisations of the BVP ([1]) that can be tackled by Trefftz methods are; 

• Neumann boundary conditions du/dn — g^ on Fo; 

• discontinuous and piecewise-constant wavenumber k; 

• piecewise constant and discontinuous tensor coefficient A in the more general 
Helmholtz equation —V ■ (AVu) — k^u = 0, e.g. llhTIl and ifTSl Ch. 1.5]; 

• spatially varying impedance 0 < tJ G LF{Fr); 

• absorbing media k gC; 

• inhomogeneous Helmholtz equation —Au — k^u = /, where the source term / 
might be either localised IIJTI §5], 024ll57ll58l . or not fT] §2.2]; 

• scattering in unbounded domains; 

• scattering by periodic diffraction gratings in II211I1 1911 ; 

• scattering by screens (i.e. manifolds with boundary, leading to non-Lipschitz 
computational domains) in 112011 . 

The presence of smoothly varying coefficients is more challenging for Trefftz meth¬ 
ods, as in general no Trefftz functions in analytical form are available; this extension 
is briefly addressed in 113.41 


in 

on fb, 


onFR. 


( 1 ) 


1.2 Notation 


We introduce a finite element partition = {^} of FI, not necessarily conform¬ 
ing. We write n^: for the outward-pointing unit normal vector on dK, and h for 
the mesh width of i.e. h := max/^g^/i^-, with Hr := diam^f. We denote by 
•^h '■= Ua'g 55, ^od := \ dFl the skeleton of the mesh and its inner part. 

We also introduce some standard DG notation. Given two elements K\,K 2 G 
a piecewise-smooth function v and vector field T on we define on dK\ n dK 2 
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the averages: {{v}} := {{t}} := 

the normal jumps: [[vJa^ := [[tJa-:= T|;fj • n/fj + . 

We denote by V/, the element-wise application of the gradient V, and write = 
n ■ V/, on do. and d„i^ = n^- • V/, on dK for the normal derivatives. 

For i > 0, define the broken Sobolev space and the Trejftz space T{£^i,): 

■= {v e : v\K e H\K) 'iK e 

T{^h) ■= {v ■ —Av — k^v = 0 in/T and dn/^v S l}{dK) \/K & ,'^h]- 

The discrete Trefftz space Vp(^) is a finite-dimensional subspace of T{^f,) and 
can be represented as Vp{^h) = 0A:e.55, ^pk{^)^ where Vpf.{K) is a pjf-dimensional 
subspace of T{^7h) of functions supported in K. We use the terms /r-convergence to 
mean the convergence of a sequence of numerical solutions to u when the mesh 
is refined, i.e. h —>■ 0, p-convergence to designate the convergence when the local 
spaces are enriched, i.e. p := min/^g^^ps: —S’ and /ip-convergence to mean the 
convergence for a suitable combination of the two refinement strategies. We remark 
that when non-polynomial spaces are used, as it is the case for Trefftz methods in 
frequency domain, it is not obvious how to define the “degree” of a space, thus pK 
denotes the local number of degrees of freedom. Finally, we denote by Re{-}, Im{ } 
and “ the real part, the imaginary part and the conjugate of a complex value. 

We note that some of the methods in 0 such as the TDG, the UWVF and the 
VTCR, involve sesquilinear forms (i.e. test functions are conjugated) while others, 
such as the DEM and the WBM, involve bilinear forms (test functions are not con¬ 
jugated). Any method (if no unbounded elements are used) can be modified to either 
form, even though sesquilinear forms are more amenable to stability and error anal¬ 
ysis; for each method we follow the conventions of the references we cite. 


1.3 Estimation of the L^(^2) norm of (piecewise) Trefftz functions 

Given two uniformly positive functions A G UTd) and cj G UTr), we 

introduce the following skeleton seminorm (defined e.g. on e > 0): 

IIMIli.a •= ll<^[[^/A]]A||i2(^/) + l|A[[v]]Ar||^2(_^/) (2) 

+ l|0’(^nV + iAt>v)||^2(f^) + ||■^v|lL2(p^) . 

A special property of the Trefftz space T{tl7h) is that this seminorm is actually a 
norm for it, and that it controls the L?{Q) norm, as it was first proved by R Monk 
and D.Q. Wang using a special duality technique in ll^ Th. 3.1]. 

Lemma 1. 111 • 11 U.o- ^ norm in Moreover, all Trejftz functions v G T {fffj n 

£ > 0, satisfy the estimate 
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IMlL2{i2) <C*|||v|||;t,a> 

with a constant C* > 0 depending only on and Setting 

Ok ;= essmf^^3K\rD := ■^(x) V/T e 

we can express the dependence of C* on the relevant parameters in the following 
situations: 

(i) If do. = Fr and Q is either convex or smooth and star-shaped with respect to a 
ball, then 


where Ci > 0 depends on d, the shape-regularity of the mesh and the shape ofQ. 
(ii)Ifk 1> 1, £2 CZ has diameter — 1 and satisfies 

x-n>7>0 a.e. on Fr and xn<0 a.e.onFo, (3) 

and each element K is star-shaped with respect to a ball of radius PRhR, we have 

1 j'2. 

+ + |||,|||,„, 

where D <it<SQ<\ 12, sq being the “elliptic regularity parameter” of 4531 
eq. (6)], and C 2 > 0 depends only on £1, i^, t, and on the shape-regularity 
infA'G^ Pk of the mesh. 

The bound in part (i) of Lemma [T] can be verified following the proof of lf85l 
Lemma 4.3.7], while that in part (ii) requires also the stability and trace estimates 
of II 54 I eq. (7), (20)] (see also ll54l Lemma 4.5] and a weaker but more general bound 
in ll5?l Lemma 4.4]). Conditions Q on the shape of £2 are satisfied if Fr is bound¬ 
ary of a domain star-shaped with respect to a ball centred at 0 and Fo is boundary 
of a smaller domain (a scatterer, or a “hole” in f2) star-shaped with respect to 0, 
see ll5^ §2, Fig. 2]. The value of the bounding constants arise only from (a) trace 
estimates for mesh elements, and (b) stability bounds for an inhomogeneous Helm¬ 
holtz BVP on £2, thus more general shapes of £2 give different dependencies on k 
(using e.g. the k-explicit//'(f2) bounds in lf30l Th. 2.4], 01061 Th. 1.6], and bounds 
in higher-order norms as in 0411 Lemma 2.12]). This result is relevant because, for 
Trefftz methods that allow a priori stability or error estimates, these are typically in 
a skeleton norm similar to ||| • |||;^, <j. Thus Lemma[T]can lead to error estimates in 
the mesh- and parameter-independent (i2) norm; we pursue this in il2.1l 112.2.11 
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2 Trefftz variational formulations 
2.1 Least squares (LS) methods 

Least squares methods are perhaps the simplest kind of Trefftz formulations. They 
allow simple error and stability analysis, are easy to implement, lead to sign-definite 
Hermitian (or symmetric) linear systems, at the price of a possibly worse condition¬ 
ing. A description of Trefftz LS schemes for the Helmholtz equation with numerous 
references is given by M. Stojek in 01071 . The same method is named frame less 
Trefftz elements in 0991 §3.6] and weighted variational formulation (WVF) in 0590 . 
In 0891 . Monk and Wang proposed the following Trefftz LS method for the BVP ([T]i: 

hnd Mls = argmin J{vhp-,gR,go), where 

'’hp^Vp{3^h) 

J{v,gR,gD)-= / fA^|[[v]]iv|^-f (7^|[[V/,v]]|^) d5 (4) 

+ / (T^|5nV + ikt?v-gffPd5-f / A^|v-gDpd5, 

Jrp ' Jro ' ' 

where [[Vvj] := on dKi n dK 2 is the jump of the complete gradient 

(whose “sign” depends on a choice of the ordering of the elements in ,^/,). The 
LS methods in 01071 eq. (7)] and OTSI Ch. 10] differ from (|4]l (apart from the use 
of different boundary conditions) in that only the normal component of the jump 
of the gradient [[V/,v]]a^ is penalised on as opposed to the entire jump [[V/,v]]. 
Obviously, every Galerkin discretisation of the variational problem arising from (|4|i 
will give rise to a Hermitian linear system, which is a clear advantage of LS methods. 

The choice of the relative weights 0 < A, <7 G between the terms in (|4]i 

is a crucial point for the conditioning and the accuracy of LS methods. Different 
choices have been proposed (for 2D problems); (7=1 and A = k or A|g = l/hg 
in §2]; A = 1 and (7|e = he/{pKi + PK 2 ) in II1071 §3.2]; A = 1 and (7|g = 
ff{max{pKi,PK 2 Y’^/^) in ITTSI Th. 10.3.4]. Here, e = dK\ PidK2 denotes a mesh in¬ 
terface, he its length, pK^ and the dimensions of the local Trefftz spaces Vp^,^ (/Ti) 
and Vpf.^{K2) on the adjacent elements K\ and K2. In 2D and 3D, li59l suggests to 
choose (7=1 and A = k and, for BVPs with singular solutions, <7|f^ = k^/^. 

The LS method computes the element Mls in Vp{^h) that minimises the error 
u — Mls measured in the skeleton norm ||v||ls := 7(v;0,0), thus orders of converge 
in this norm follow immediately from approximation bounds for the specific discrete 
Trefftz space Vp{£Tp) chosen, see e.g. Sj^below or ll89l . Since 11 |v| | |;l ^ <I|v|Ils (with 
equality if 7 in (|4|i is defined with [[V/,v]]aj instead of [[V/,v]]), Lemma[T] following 
lf89] Th. 3.1], guarantees that the L?'{Q) norm of the error of the LS solution is 
controlled by the value of the LS functional, thus convergence follows also in Q. 
This is summarised in Theorem[T] see ill.3l for the extension to different domains. 
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Theorem 1. Let u be the solution of ([TJ and ^LS ^ Vp{^h) the discrete LS solution 
of®. Then, for C* > 0 depending only on k,X,G,Td,Q and 






|m-v/,p||ls- 


IfX = k, G = I, dQ = Fr and Q is either convex or smooth and star-shaped, then 


! —MLs||z,2(i2) ^ Co diami 2 A: ^/^(^ + (kmmhK) 


inf M —V/ 2 P 

Vhp^^pi^h) 


where Co > 0 depends only on d, the shape of Q. and the shape-regularity of 
The /ip-convergence theory of 1541 easily extends to the LS method. In 2D, if the LS 
parameters are defined as Xj^ = kh/m\n{hK^^hK2} for o = dK\ n dK2, X^^ = kh/hR 
for e C dK flLb, and G^ = l/k, under the assumptions on Q. and on the discretisation 
stipulated in ll54l . then the || • || norm of the LS error is estimated as in ll54l eq. (48)] 
and the Lf{Q.) norm of the same error converges to zero exponentially in the square 
root of the total number of degrees of freedom used. 

In ll75l Ch. 10], the Trefftz LS scheme is analysed for pure Dirichlet boundary 
conditions (Lr = 0 ); the crucial parameter in the analysis is the relative distance 
between k^ and the closest Dirichlet eigenvalue of —A. Error bounds in the broken 
Sobolev norm are derived. 

In the numerical tests in and li40l . the LS method appears to be slightly less 
accurate than the UWVF (see il2.2.2l below') and a DG method, all employed with 
the same discrete space. On the other hand, in the examples in li59]| , the performance 
of the LS method is comparable to that of the UWVF and considerably better than 
that of the VTCR. 


2.1.1 The method of fundamental solutions (MFS) 

A popular class of LS Trefftz methods is the method of fundamental solutions. A 
lucid introduction to the MFS for Helmholtz problems, together with numerous ref¬ 
erences, is in lUl . The MFS is considered a special case of source simulation tech¬ 
nique in |l92l. The characteristic features of the most common form of the MFS are: 
(i) the domain is not meshed; (ii) the N basis functions are fundamental solutions 
(//g*^(k|x — y^l) in 2D, i = where Hq'^ is a Hankel function of the first 

kind and order zero and G \ Q, see il3.3l) : (Hi) the minimisation of the Lf{d£2) 
norm of the error is substituted by the minimisation of the squared error overM > N 
points Xj G dQ, j = 1,... ,M. If M = V, the MFS is not an LS method but it simply 
interpolates the boundary conditions with Trefftz functions. 

The same method with plane wave bases is compared to the MFS in H]. A variant 
that is popular in acoustics is the Helmholtz equation least-squares (HELS) method, 
which uses spherical-wave and multipole basis functions, see the recent book 111 171 
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and references therein. LS variants of MFS relying on higher order multipoles in 
addition to simple Hankel functions have a long history in wave simulations ll90l §2]. 

The locations of the basis singularities are either obtained numerically together 
with the coefficients multiplying the basis functions using non-linear LS solvers lISTl 
eq. (7)] (leading to a highly adaptive method), or can be fixed a priori on a smooth 
boundary in M" \ Q, e.g. using complex analysis techniques (in 2D) as in 19], or are 
determined based on heuristic criteria 1^ §3]. 

The MFS with fixed nodes can be interpreted as a discretisation of a compact 
transfer operator related to a single layer potential representation. For this rea¬ 
son it yields ill-conditioned linear systems; however this does not rule out effi¬ 
cient computations as demonstrated and analysed in @ and in ifTOl §7]. Accord¬ 
ing to Il3n p. 766], the larger the distance between the nodes and Q, the more ill- 
conditioned the linear system and the more accurate the solution (though this might 
seem counter-intuitive). 

A strength of the MFS is its simplicity of implementation, as no mesh is needed 
and all geometric information is contained in only two point sets C K" \ f 2 , 

{^j}f=i C dQ. Since fundamental solutions satisfy Sommerfeld radiation condi¬ 
tion, the MFS is often used for scattering problems in unbounded domains. 

In 121, the convergence of the MFS for Dirichlet problems on a circular domain 
is analysed in great detail, and a special design of the curve supporting the funda¬ 
mental solutions is proposed for general domains with analytic boundaries. With 
this choice, extremely accurate and cheap computations are possible. 

In ITOl . Barnett and Betcke present a finite element scheme that couples the LS 
formulation of II 1071 with the MFS in 2D. They consider the scattering by sound-soft 
(non-convex) polygons; the total field is approximated inside an artificial boundary 
and the scattered field outside of it. Singular Fourier-Bessel basis functions de¬ 
pending on the scatterer’s corners (see 0.41 i are used on all elements adjacent to 
the scatterer, strongly enforcing the (homogeneous) Dirichlet boundary conditions; 
due to this, no terms on dQ appear in the method formulation. Exponential orders 
of convergence are proved. The strong enforcement of boundary conditions may be 
substituted by an LS approach to deal with more general linear boundary conditions, 
curved boundaries and transmission problems. 


2.2 Discontinuous Galerkin (DG) methods 

The discontinuous Galerkin (DG) methods constitute a wide class of numerical 
schemes for the approximation of PDEs, employing discontinuous test and trial 
functions i). A great number of tools for their design, implementation and error 
analysis have been devised, so they are a natural setting for Trefftz methods. In ll55l 
we showed that when the interior penalty (IP) method, one of most common DG 
schemes, is applied to the Laplace equation, the use of Trefftz spaces (made of 
harmonic polynomials) offers better accuracy than standard spaces also in an hp- 
context. Similar considerations were made in lEl for the /i-convergence of the local 
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DG (LDG) method. To our knowledge, no standard DG variational formulation (e.g. 
any of those in |0) has been proposed in the literature to discretise time-harmonic 
problems with Trefftz basis functions. Possible reasons for this are that the error 
analysis of standard DG schemes requires inverse estimates, which are well-known 
for polynomial spaces but harder in the Trefftz case (however, see ll46l §3.2] for 
/t-explicit inverse estimates for plane waves in 2D), and that the application of for¬ 
mulations designed for the Laplace equation to the Helmholtz case requires some 
problematic minimal resolution condition to ensure unique solvability ll82l . 

In the next subsections we outline some DG formulations that have been designed 
specifically for Trefftz discretisations; some of these have later been employed also 
with polynomial approximating spaces, e.g. 082II88I . 

A note on terminology: all Trefftz methods presented in this survey involve the 
discretisation of variational formulations based on discontinuous functions, how¬ 
ever with “DG” we denote only those methods that arrive at local variational for¬ 
mulations by applying integration by parts to the PDE to be approximated. On the 
contrary, least squares and weighted residual methods simply enforce (weakly) con¬ 
tinuity and boundary conditions, irrespectively of the considered PDE. 


2.2.1 The Trefftz-DG (TDG) method 


Originally, Trefftz-discontinuous Galerkin (TDG) methods (or plane wave DG, 
PWDG, when used in combination with plane wave basis functions) were in¬ 
troduced as a way of recasting the ultra weak variational formulation (UWVE) 
of lfT8][T9l (see il2.2.2l below) in a framework that would facilitate its theoretical 
analysis inillll. A similar, but more general, Trefftz-DG framework was proposed 
in II37II39I . arising from methods for hyperbolic equations; see Remark[T]below. 

We first derive the TDG formulation as in ll5^ . We multiply the Helmholtz equa¬ 
tion O by a test function v and integrate by parts twice on each K & 

0 = [ {—Au — k^u)vdV = f (Vm • Vv —k^Mv)dy — [ Vu nKvdS 
JK Jk JdK 

= / u{—Av — k^v)dV+ / udn^vdS— / d„j.uvdS. 

Jk JdK JdK 

We then replace u and v by discrete functions u^p^Vhp €Vp{3J'k), the trace of u on 
dK by the numerical flux M/,p, and the trace of Vm by the numerical flux ikGhp (both 
defined below), obtaining the elemental TDG formulation: 


/ 

JdK 


^hp 




I 

JdK 


ikoi^p • v/jp diS — 0, 


(5) 


where the volume integral vanishes as the test function v/,p GVp{J^h) C is a 

Trefftz function. Variants of DG methods are distinguished by the underlying nu¬ 
merical fluxes. Here we opt for the primal fluxes: 
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{ {{V/jW/jp}} (X\k [[w//p]]/i/ 

^h^^hp (1 h^hp “I” ik'O^UppH 

^h^hp CCik {lij^p 

f {{«/.p}} - P V^hUhph 

Uhp= Uhp-5{{ik-&y''^hUhp-n + Uhp-{ik-&y'gR) 
[gD 


on faces in 
on faces inf];, (6) 
on faces in fb, 

on faces in 
on faces in Fr, (7) 
on faces in Fd, 


where the flux parameters a > 0, j3 > 0, 0 < 5 < 1 /2, are bounded functions defined 
on suitable unions of edges/faces (see also Table [T]l. Adding over all elements, we 
obtain the following formulation of the TDG method; 


And Mtdg € Vp(,^,) s-t. = 2 /tdg(mtdg, v/,p) = £tdg(v/,p) Vv/,p G yp(5ft), where 

■*2^dg(m,v’) := (8) 


({{«}}- {{^/2 «}} • Ma' + «i^[[M]]Af • [[vjiv - P (ik) ^ huhl'^hv]]N^ 

+ J ^(1 — 5)i^t?Mv+(1 — 5)M(9nt'— V—5(iA:jJ)^^(9nM(9nV^ d5 
+ J ^ — (9nM v+cti^Md5, 

^tdg(v) := gR(^{l - 5 )v- 5{ikT^y^d„v^ Jr SD(^a\kv-dn 


d5 


' — <5nV d5. 


The TDG method was introduced in the primal form described here in II44I46I and in 
mixed form in lf56l . under the name of plane wave DG (PWDG) method, following 
the derivation of lb] of general DG schemes for elliptic equations. In ll4^ . first- 
order convergence in the meshwidth was established, using Schatz’ argument, for 
2D Robin problems with source term / G L^{D,), plane wave discrete spaces and 
quasi-uniform families of meshes. This was extended to higher orders in h in ll84l . 
p-convergence in lf52l . three dimensions in ll85l . locally-refined meshes in ll5^ . and 
Anally the exponential convergence in the number of degrees of freedom of its hp- 
version was proved in ll54l . Its dispersion analysis was performed in 04411451 . 

For polynomial discrete spaces, the advantages of using the formulation under¬ 
lying the TDG method, compared to standard DG schemes, were analysed in ||82l. 
In ifTSl . the TDG formulation was utilised with (non-Trefftz) bases defined from os¬ 
cillating functions from high-frequency asymptotics modulated with polynomials; 
problems with varying coefficients were also considered. 

The TDG formulation ® can be seen as a modification of either the interior 
penalty method, or of the local DG (LDG) method (see e.g. ||6l): with respect to 
the interior penalty method, the stabilisation term multiplied by j3 is added in the 
TDG fluxes 0, while with respect to the LDG method, in the TDG fluxes (|6]l, the 
consistency term is written in terms of the primal variable ({{VpMpp}}) instead of 
in terms of the auxiliary variable ({{i^CJ/ip}}) and the additional stabilisation of the 
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jumps of Ghp is removed. In Il05l . the TDG and the UWVF are seen as special 
instances of a family of methods arising from integration by parts. 

The a priori error analysis of the TDG relies on Theorem|2]below (e.g. ||5^ §4]), 
which makes use of the following mesh- and flux-dependent seminorms; 


:= k 




+ k-^ 

I |||2 _ III |||1 

r 11 Itdg+ ■ 11 r 11 Itdg 

+ k 


^dnV 


lHH) 
2 


+ k 




lHH) 




H” k 


(1 




+ k 




2 


+ k 




P-HM} 

2 




+ k-^ 


LHrp) 




a ^dn\ 


a-3{{V,v}} 


lHFd) 


lHH) 


Theorem 2. The seminorms 111 -11 |tdg o.nd 111 • 11 |tdg+ norms in the Trefftz space 
T{ffff). The TDG sesquilinearform is continuous and coercive: 

Ktdg(v,w)| <2|||y|||T.DG+llk|||TDG, Im{ jz/tdg(V, v)} = |||v|||tdg 

for all V, w G T{^i,), thus there exists a unique solution Mtdg S Vpif^h) to the TDG 
formulation i} and the quasi-optimality bound holds: 

11 |m — MtdgI I Itdg ^ 3 inf 11 |m ~ r/,„| | |^j3q+. 

Vkp^VpiSu) 


Choosing = ak on UTb, = P/k on and <7^ = min{5,l — 5}/2kt> 
on Fr, the norm @ is controlled as |||v|||;l < |||v|||tdg for all v S T{tffji)- Thus, 

by Lemma[T] the X{£2) norm of the TDG error can be controlled by its ||| • |||tdg 
norm, and so by the discrete space approximation properties. This result has been 
stated in several slightly different forms, depending on the regularity of the solu¬ 
tion M, the type of mesh used, the choice of the numerical flux parameters a,(5,5', 
see ll^ Lemma 4.3.7], ll5^ Lemma 4.4] and ll54l Lemma 4.5]. To strike a balance 
between the size of the constants arising from the duality argument of Lemma [T] 
and approximation errors, different flux parameters have been chosen on different 
meshes and aiming at different types of convergence estimates, see Table [1] For il¬ 
lustration, we state the result in the case of constant flux parameters, quasi-uniform 
meshes, and domains that guarantee sufficiently smooth solutions for the dual prob¬ 
lems; this follows from Lemma[T]and Theorem|2](c/; 185] Cor. 4.3.8]). 

Corollary 1. Let u be the solution of ([T]l, where Q. is either convex or smooth and 
star-shaped, and let Mtdg C be the solution of the TDG method with flux 

parameters chosen as in the second row ofTable\I\ Then 

||M-MTDG|lz, 2 (i 2 ) <Codiamf2 (l + (kmin/r/j:)^^/^) inf |||m-v/,p|||tdg+, 
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where Cq > 0 depends only on d-, the shape of Q and the shape-regularity of the 
mesh, but is independent ofk and Vp(^). 

The combination of the abstract error analysis outlined above and approximation 
estimates for plane, circular and spherical waves (see leads to a priori h-, p- 
and /ip-convergence estimates in ||| • |||tdg and if norms, see 046115245^[85]| . The 
dependence of the error bounds on the wavenumber k is explicit, as in Corollary[T] 



a 

P 

S 

Quasi-uniform meshes, Ij-convergence 

|46] 

a.jkhK 

hkliK 

dkhx 

Quasi-uniform meshes, p-convergence 

(52] 

a 

b 

d 

UWVF fsee 

dH 

1/2 

1/2 

1/2 

Locally refined meshes, /ip-convergence 

ED 

ah/hx 

hh/hx 

d/j / hx 

Geometrically graded meshes, exponential Itp-convergence 1541 

ah/iiK 

b 

d 

Polynomial (non Trefftz) basis, hp-convergence 

(82) 

aql/khx 

bkhx/qx 

dkhx/qx 


Table 1; Different TDG flux parameters in (| 6 ]l and (|7]i that have been considered. 
Here a,b,d are positive functions independent of the other parameters; k is the 
wavenumber; hK is the local meshwidth; h — maxjfg^j^ hx is the global meshwidth; 
qK is the local polynomial degree (for the non-Trefftz version). 


Remark 1. The Helmholtz equation may be written as the first order hyperbolic 
system —i^u + = 0 , where u := {u^Vu/{\k)) and Af-') are the 

(1 +n) X (1 +«) symmetric matrices whose only non-zero elements are = 

1 = 1, for 1 < j < n. Then, similarly to llJTl eq. (22)] or ||^ eq. (5)], a general 
Trefftz-DG method can be written as: 

seek u e Vp(,3^,) := {(m,(t) : m S Vp(,5^,),(J = VM/(ik)} s.t. VveVp(5^,) 

L L _ (F“u-g)-vd5 = 0 

Ki^K2 

where the flux-splitting matrices F™, F°“'^ are defined on YIkgSi, and satisfy F™ < 
0 , F™‘ > 0 (i.e. are negative and positive semi-definite, respectively), F™ -f F°“* = 
() on dK, and F^^ = —on dKi n dK2. The boundary data are represented 

by a suitable vector held g = —F°“‘u. The TDG in (HJ (up to a factor —ik) is obtained 
by choosing: 
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-a 


inT 


in/f -pn^n' J 
—(1 —5)jJ 5nJ 
(l-5)n -l-iKg)] 

-«nA 

0 0 


F ont _ 

K - 



on dKn^l, 
on dKCiFu, 
on dKCiFo. 


The right-hand side is represented by the vector g = on and 

g = -(n“ kDonTo. 


2.2.2 The ultra weak variational formulation (UWVF) 


The ultra weak variational formulation (UWVF) has been introduced in the 1990’s 
by O. Cessenat and B. Despres in II18II19II . Since then it has received a great deal 
of attention and has been applied to numerous PDFs and BVPs; we refer to ED 
for a description of its computational aspects and to ff^ §3.5.2] for an extensive 
bibliography. Different derivations can be found e.g. in 0 1 711 1 9|[J7][^l4^ : in partic¬ 
ular HI711461 obtain the UWVF in the setting of DG schemes for elliptic problems 
of El, while IIJTII^ derive it for general first-order hyperbolic systems using a 
flux-splitting approach as we did for the TDG in Remark[D Note that different pa¬ 
pers use different sign conventions. The extension of the UWVF to problems with 
smooth coefficients has been tackled in ll65l . 

To write its formulation for the BVP ([D in the Robin case, i.e. Tb = 0, we first de¬ 
fine the trace space X L^{dK), and the operators Fk : L?{dK) L^{dK), 

mapping the boundary datum yx of a local adjoint-impedance Helmholtz BVP into 
the impedance trace of the BVP solution ck itself: 


FK{yK) ■■= {dtiK + ik)eK, where 


—AeK — k^eK = 0 in/f, 
{-dtiK + ^k)eK = yK on dK. 


The Helmholtz BVP is written as a transmission problem across the mesh interfaces, 
i.e., for all K,K' e ,5/,, 


—Au —ku = 0 in/f, 

d„f,u + iku =—d„^,u + iku ondKCidK', 
d„j^u + ik'&u = gR ondKnFR. 

Then, after multiplying the first equation by e G T{£^h), integrating by parts 
twice, taking into account transmission and boundary conditions, and introducing 
x,y GX defined as = —dn^^u + iku and y^^R = -f ike, the UWVF of prob¬ 

lem ([D lfT9] (1.4)] reads: find xGX such that, for every y GX, 
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dKnTR 1 + J? 


2 


■8RFK{y\dK)dS. 


(9) 


(Note that for t? = 1 the term on dKnPR at left-hand side vanishes and 2/(1 -f t?) = 
1.) The expression (|9l) is a variational formulation for the skeleton unknown x; after 
the equation is solved for x, the Helmholtz solution M|^ can be recovered in the 
interior of each element by solving a local (in K) adjoint-impedance Helmholtz 
BVP with datum -f ik)u^K = x\;)r. If the formulation is discretised choosing 

a hnite dimensional subspace X/, of X corresponding to the impedance traces of a 
Trefftz space, namely 


Xji {^Xh S X ; Xh\^gR — {—dn^ -\-ik)v\R VK £ v £ Vp{3^h) } 


then the action of Fr and the reconstruction of ur in K are immediately computed. 

Theorem 2.1 of ||T91 states that the discrete problem obtained by substituting X/, 
to in (|9l) is solvable, independently of the meshsize /i; Corollary 3.8 shows that, 
for plane wave discrete spaces, the Dirichlet and Robin traces of the UWVF solution 
converge to the corresponding traces of u with algebraic orders of convergence in 
L^{rR). In IITtI §4], these results have been used together with the duality technique 
of Il89l to prove orders of convergence for the L^{Q.) norm of the error. 

The UWVF has been recast as a DG method with Trefftz basis functions in 
several different ways in In particular, 1461 Remark 2.1] shows 

that the UWVF is a special case of the TDG formulation ® for flux parameters 
a = P — 5 — 1 /2. As a consequence, the orders of convergence in h and p proved 
for the TDG on quasi-uniform meshes in l46lf52l carry over to the UWVF (with 
suboptimal orders in h); on the other hand, the hp-type results of l5^l54l require 
variable numerical flux parameters to cope with elements of different sizes (see Ta¬ 
ble [T]i, so they do not apply to the UWVF. Thus, the TDG can be understood as the 
extension of the UWVF to non quasi-uniform meshes. Alternatively, in l88l §4.3, 
5.2], the UWVF is employed on meshes refined towards solution singularities by 
choosing Trefftz spaces on large elements and polynomial spaces on small ones. No 
applications of the TDG combining mesh-dependent parameters and polynomial 
spaces in small elements have been documented. 

2.2.3 DG schemes with Lagrange multipliers 

The DG schemes described so far enforce weak continuity between elements using 
numerical fluxes, in the spirit of |0. A different approach is to enforce continuity 
using Lagrange multipliers. This was probably first proposed for Trefftz methods 
in ll63l §2.3], for the ID Helmholtz equation. 

This strategy has been followed in the discontinuous enrichment method (DEM), 
introduced by C. Farhat, I. Harari and L.R Franca in ]|32|, combining a space of 
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piecewise-constant Lagrange multipliers on mesh interfaces with a discrete space 
composed by sums of continuous piecewise polynomials and discontinuous plane 
waves. Subsequently, in ll^ . the polynomial part of the trial space was dropped, 
leaving a plane wave trial space and thus reducing to a Trefftz method; in this ver¬ 
sion, the DEM was renamed discontinuous Galerkin method (DGM) and the La¬ 
grange multipliers were approximated by oscillatory functions. This formulation 
performed very well for test cases and was later extended to “higher order ele¬ 
ments” (i.e. elements containing more plane waves) and other PDEs. We refer again 
to ll76l §3.5.3] for a comprehensive bibliography. 

Here we briefly describe the formulation of the DGM following ll3?l §2]: 


find (m,A) G xW{3ih) s.t. 



where 




It is immediate to verify that the solution u to BVP ([T]i satisfies this formulation, and 
that the multiplier A represents the normal derivative of u on the mesh interfaces 
and on Tq. This formulation is then discretised by restricting it to finite dimensional 
spaces Vp{£^h) C and Wp{£^h) <zW{^h). In the DEM of ([321, Vp{^h) is the 

direct sum of a continuous polynomial and a plane wave space, in the DGM of 1^ 
and subsequent papers only the plane wave part is retained, soVp{^h) GT{£^i,). The 
volume degrees of freedom, i.e. those corresponding to Vp{^j), are then eliminated 
by static condensation in order to reduce the computational cost of the scheme. 

A stability and convergence analysis of the simplest version of the DGM (four 
plane waves per element and piecewise-constant multipliers) is attempted in 0: 
for a Robin-Neumann BVP on a domain decomposed in rectangles, under a mesh 
resolution condition, the scheme is shown to be well-posed, and a priori orders of 
convergence are proved (in //' (,5^,) norm for the primal variable and in for 

the multipliers), along with residual-type a posteriori error bounds. We are not aware 
of any error analysis for the DGM method holding in more general situations (e.g. 
more than four plane waves per elements, propagation directions not aligned to the 
mesh, non-rectangular mesh elements). 

A similar formulation, named hybrid-Trejftz finite element method, is described 
in ll99] §3.5] (deriving the functional in eq. (65) therein); the same form £/dgm 
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above is used, while SSoam is substituted by := — fji-/ H [[V/,w]]ivd5 — 

II 5nwd5, where now the multiplier jj. approximates the Dirichlet trace of u, the 
right-hand sides and the space W{^h) changed accordingly. A further variant of 
hybrid-Trefftz methods is presented in Ml091 and related papers. 

Another DG method with Trefftz basis, called modified DG method (mDGM), 
has been proposed in ll48l . The Lagrange multipliers are double-valued on the in¬ 
terfaces (differently from the DEM/DGM of M32II33I 1 and belong to \ 

Fr). a two-step procedure is adopted. First, for each basis element A € L?{dK\rR) 
of the discrete Lagrange multiplier space, a well-posed Helmholtz BVP on K with 
impedance datum A is solved in the local Trefftz space Vpp{K) using the classical 
(/T)-conforming variational formulation. Second, these local solutions are com¬ 
bined in a global LS formulation leading to a positive semi-definite system whose 
unknowns are the Lagrange multipliers themselves. The mDGM was further im¬ 
proved in ||2l leading to the stable DG method (SDGM), which differs from the 
mDGM in that the local impedance problems are solved with a least squares formu¬ 
lation posed on dK, which gives local Hermitian matrices. 

Lagrange multipliers are also used to tackle problems with discontinuous coeffi¬ 
cients by means of the partition of unity method, see II 73 I and il2.5l below. 


2.3 Weighted residual methods 

Trefftz discretisations lend themselves well to weighted residual formulations: the 
discrete solution is automatically a local solution of the PDF, only the residual on 
interfaces (the jumps) and on the boundary (the mismatch with boundary conditions) 
need to be enforced by multiplying them to suitable traces of test functions. The 
choice of these traces leads to different variational formulations, the most developed 
of which are the VTCR and the WBM described in the following. While it is simple 
to design weighted residual methods, their error analysis is by no means easy, as 
they arise neither from integration by parts, nor from a minimisation principle. 

An earlier weighted-residual Trefftz formulation is the weak element method of 
ill, where the integral averages of Dirichlet and Neumann jumps on mesh faces 
are set to zero (equivalently, test functions are constant on each mesh face). 

We note that some of the earliest Trefftz schemes, e.g. the indirect approximation 
of fSA eq. (35)], are of weighted-residual type, even though testing was confined to 
the boundary of the domain only, see il2.4l below. 


2.3.1 The variational theory of complex rays (VTCR) 

The VTCR is a weighted residual Trefftz method introduced in the 1990’s by 
P. Ladeveze and coworkers for problems arising in computational mechanics and 
later extended to the Helmholtz case in fllOll . Recent surveys are M70ll71IIT00l . 
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Several VTCR formulations, slightly different from each other, have been pre¬ 
sented. A general VTCR formulation for the BVP ([T]i can be written as: 


find Mvtcr C Vp(,5/,) s • t. .2^TCr(MvTCRj t'/ip) — ^VTCR {vhp) '^Vhp e Vp(^,): where 


M'tcr(m,v) :=Im 


(NW ■ {{VW} - MivM}) 


d5 


( 10 ) 


+ [ M(9nvd5-|- f (■^^{daU + ik'd-u)dj,v + C2{daU + ik'd-u)v\ dS 
Jrn Jr,, ViKj> / 


^vtcr(v) Imj^ gDdnVdS + (^■^gRdaV + C2gRV^dS 


where we have reported the formulation with only the imaginary part of the left- and 
right-hand side, following the VTCR convention; however dropping ’Tm” does not 
modify the method. 

The formulations in 01001 eq. (21)] and in ll?!] eq. (5)] correspond to the choice 
of coupling parameters Ci = 1/2 and C 2 = —1/2 (up to an overall factor k and 
using Re{—/z} = Im{z}); that in 01021 eq. (6)] to Ci = 1/2 and C 2 = 1/2; that 
in ll68] eq. (4)] to Ci = 1 and C 2 = 0. The choice of the coupling parameters does 
not affect the consistency of the method as all terms in (fTOl l are products of residuals 
(internal jumps and boundary conditions) and traces of test functions. In some of 
the papers cited, using lm{ab} = — lm{ab}\/a,b € C, the conjugation is written 
on the trial, rather than test, functions in some of the terms, without affecting the 
formulation. 

The VTCR (and similarly the WBM) does not correspond to any of the classical 
DG schemes listed in Jb). Indeed, to derive it from the elemental DG equation (|5]l, 
one would need to choose numerical fluxes that, in the terminology of lb), are nei¬ 
ther consistent (they do not equal the helds Vm and u when applied to the exact 
solution u itself) nor conservative (they are not single-valued on the interfaces). 

Following ll68] §2.2], it is possible to show that if absorption is present then the 
VTCR is well-posed. More precisely, provided that Ci = 1, C 2 = 0, ReA: > 0 and 
Im{A:^} > 0, the VTCR bilinear form satishes 


■*2^TCr(v’, v) 


-Im{A:^}||v||^2(i2) 


Re;t 

w 




2 

LHrp) 


'^veTi^h), 


thus the VTCR solution is unique in the Trefftz space and coercivity in L^(X2) norm 
holds (the analogous result for Ci = —C 2 = 1 /2 is proved in 11711 Prop. 2]). However, 
this does not extend to the setting we considered so far, i.e. propagating waves with 
A: G R: in this case it can easily be shown that .2 /vtcr('', v) = 0 for all v G such 

that V = 0 on all elements adjacent to the Robin boundary Fr and for any choice 
Cl, C 2 G C, thus well-posedness can not be ensured using a coercivity argument. 
Following E] Prop. 2], for Ci = 1/2,C 2 = -1/2,/t G M, we have: 


^VTCr(v, v) — — — t? 


C{r«) 


+ k 






VvGr(,^,), 
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thus (using Holmgren’s theorem Il20l Th. 2.4]) uniqueness of the solution of ( fTOb is 
proved if all mesh elements are adjacent to For more general cases, coercivity 
appears to be too strong an argument. We conjecture that discrete inf-sup conditions 
might be a more viable way for proving well-posedness of the VTCR. 

Section 3 of CD considers the application of the VTCR formulation, cor¬ 
rected with suitable volume terms, with non-Trefftz (piecewise-polynomial) discrete 
spaces. This variation is termed weak Trejftz and analysed therein. 

2.3.2 The wave based method (WBM) 

The WBM is a weighted residual Trefftz method, analogous to the VTCR, first 
introduced in the dissertation of W. Desmet 12^ and later extended to a wide variety 
of engineering applications, mainly in the realm of vibro-acoustics. Recent reviews 
of the state of the art of the research on the WBM can be found in ll24llZ7l . The 
discrete space typically used together with the WBM is composed of propagating 
and evanescent plane waves, as outlined in ^3.21 

The basic variational formulation of the WBM applied to BVP ([T]|, translating 
§4.1.4 of IZTl to our notation and multiplying all terms by {—ik), reads 

find Mwbm G Vp{^h) s .t. .!2 ^WBm(MwBM: Vhp) = ^WBM (vhp) Vv/,p G Vp{3h), where 



where Z,„, is an interior coupling factor. In some works, a slightly different formula¬ 
tion is used, e.g. in ||98] eq. (81)] different terms are used on the internal interfaces. 
We are not aware of any rigorous stability or error analysis of the WBM formulation. 


2.4 Single-element direct and indirect Trefftz methods 

Most schemes described so far were introduced not earlier than mid 1990’s, but 
a lot of research on Trefftz methods has been carried out since the late 1970’s by 
I. Herrera, J. Jirousek, A.R Zielinski, O.C. Zienkiewicz and numerous co-workers, 
mainly for static elasticity problems. General reviews of these works are in 067I121II ; 
the Helmholtz case is described in detail in ll22l . A major difference between these 
methods and those we described in the previous sections is that in many instances of 
the former ones no mesh is introduced on the domain Q, so that the unknowns are 
defined on dQ only. For this reason, these Trefftz methods more closely resemble 
standard boundary element methods rather than finite element schemes. 
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There are two main classes of these Trefftz methods: direct and indirect. (We use 
the terms “direct” and “indirect” as in 022II67I and 1981 §5.1].) We describe them for 
a modification of BVP ([T]i where we drop the Robin boundary Tr and we consider 
instead a Neumann boundary portion Fn with boundary condition d^u = g^. 

The indirect method is the simplest kind of weighted residual scheme: 



( 11 ) 


(see 1^ eq. (35)] for sound-hard scattering problems in unbounded domains, l98] 
eq. (47)], 01211 eq. (16)], WT\ eq. (16), (26)]). For Dirichlet exterior problems this 
is also the method of 0 §3]. In most references the test function is not conjugated. 
We note that the indirect method is nothing else than the WBM of 112.3.21 posed 
on a single element, i.e. = {£1} and = 0. In the indirect method, the trial 
functions approximating u are global solutions of the Helmholtz equation on the 
whole of on the other hand the test function v only needs to be defined on dFl. 
If the Trefftz test and trial spaces coincide, then the obtained stiffness matrix is 
symmetric (by Green’s second identity). If the signs of the terms on Tv are changed, 
as in 1 ^ eq. (22)], a non-symmetric formulation is obtained. 

Subtracting from (fTTl l the second Green’s identity Jg^iudaV — (9nMv)d5 = 0, 
which holds for all Helmholtz solutions u and v in £2, we derive the direct method: 



( 12 ) 


(see li22l eq. (42)], li9^ eq. (50)]). The direct method for the Dirichlet problem may 
be viewed as the TDG of 112.2.11 with a = 0 posed on a single element K — FI. 
Conversely to the indirect method, consistency of (fT^ is guaranteed only if the test 
functions are Helmholtz solutions in FI, while the trial functions might be defined 
(and often are) on dFl only, for better computational efficiency; the solution is then 
evaluated in FI with a representation formula in a post-processing step as for BEMs. 
The stiffness matrix arising from the direct formulation (fTSI) is the transpose to that 
of the indirect method (fTTI) . Theorem 6.44 in 11051 gives sufficient conditions for 
the well-posedness of the direct method. Theorem 7.19 in IItTII proves that, for well- 
posed Dirichlet problems with [dFl ) data, if the Neumann traces of the trial space 
coincide with the Dirichlet traces of the test space, then the direct method is well- 
posed and computes the best approximation of the exact solution in L?{dQ) norm. 
If Q is unbounded, the direct and the indirect methods can still be used choosing 
discrete functions that satisfy Sommerfeld radiation condition; however in (fT^ the 
conjugation on the test function must be dropped to preserve consistency. In this 
case, if a multipole basis is used. Waterman’s null-field method is obtained, see llTSl 
Ch. 7], which is a special instance of the T-matrix method ifTSl §7.9]. (Note that li92l 
uses the name null-field method for the indirect method with non-conjugated test 
functions, and Cremer equations for the same with conjugated test functions.) 

For a special choice of Trefftz test functions v indexed by a complex param¬ 
eter (see the last paragraph of ^3.21 i. method (fTSli is called ‘"global relation” and 
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is the variational formulation at the heart of the Fokas transform method, see f23\ 
eq. (2)], II 1051 eq. (6.142-143)] or 11211 eq. (7.156)]. In this context, this formulation 
is typically discretised using piecewise-polynomial (on dQ) trial functions, even 
though Trefftz functions may be used as well. 


2.5 Non-Trefftz methods with oscillatory basis functions 

The main reason for the success of Trefftz methods in the context of time-harmonic 
wave problems is that the oscillatory basis functions may offer much better approx¬ 
imation properties than piecewise polynomials used in standard FEMs. On the other 
hand, similar approximation can also be achieved if the discrete functions are not 
exact local solution of the PDE to be discretised, but are are only “approximate so¬ 
lutions”. If basis functions of this kind are used, the Trefftz formulations described 
in the previous sections cannot be employed as they stand, because the residual in 
the elements will not vanish any more and consistency will fail. 

Approximate Trefftz functions are especially attractive for problems with smooth¬ 
ly varying material parameters, where no analytic Trefftz function might be known. 
Trefftz formulations, possibly with additional volume terms, can be used with ba¬ 
sis functions that are solutions of the equation only up to a certain order; see 
csiisiini, where this idea is pursued for DG, UWVE and DEM formulations. 

In the following we briefly discuss a few methods that have been proposed em¬ 
ploying oscillatory and k-dependent basis functions that are not Trefftz. 

A very well-known scheme of this kind is the partition of unity method (PUM or 
PUEEM), introduced by I. Babuska and J.M. Melenkin the mid 1990’s, see e.g. lISTll . 
The PUM combines the approximation properties of Trefftz functions with the stan¬ 
dard variational formulation of the problem, e.g. for the BVP ([T]i with Fp — 0 




This requires the use of //*(f2)-conforming trial and test functions, thus continuity 
on interfaces needs to be enforced strongly, which is not viable in Trefftz spaces. 
The PUM uses as basis a set of Trefftz functions multiplied to a partition of unity 
defined on a EEM mesh, e.g. piecewise linear/multilinear polynomial EEMs on sim- 
plicial/tensor elements. Theorem 2.1 in ifSTI ensures that the trial space obtained 
enjoys the same approximation properties of the Trefftz space employed. If a p- 
dimensional local Trefftz space is used in each element, together with a piecewise 
linear/multilinear partition of unity, the total number of degrees of freedom used 
equals p times the number of mesh vertices, while for a similar Trefftz method on 
the same mesh (providing comparable accuracy) it would equal p times the num¬ 
ber of mesh elements; this means that on tensor meshes almost the same number 
of DOFs would be employed by the two methods, while on triangles and tetrahedra 
a saving of a factor up to two or six, respectively, can be achieved by the PUM. A 
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shortcoming of the PUM is that the formulation (fT3l l is not sign-dehnite and its well- 
posedness requires a scale resolution condition, while this is not needed for some 
Trefftz schemes such as the TDG/UWVF presented in il2.2.1l and ^2.2.2\ Differently 
from Trefftz schemes, the implementation of the PUM requires the computation of 
volume integrals; moreover, the numerical integration of the PUM basis functions 
may be more expensive than that of genuine Trefftz functions, see il4.1l 

The PUM for the Helmholtz and other frequency-domain equations was further 
developed by RJ. Astley, P Bettes, A. El Kacimi, O. Laghrouche, M.S. Mohamed, 
E. Perrey-Debain, J. Trevelyan and collaborators, see e.g. 072II96I . When a PUM 
and a standard EEM discrete spaces are combined, e.g. using formulation (fOl l. the 
method obtained is termed generalised finite element method (GEEM); e.g. Coll 
employs high-order tensor-product polynomials summed to products of plane waves 
and bilinear functions. In problems with discontinuous wavenumber k, the PUM can 
be applied by coupling the homogeneous regions by means of Lagrange multipliers 
as in GS; this is not necessary as formulation (fOT l holds on the whole domain, but 
enhance the accuracy as in each subdomain only basis functions oscillating with the 
correct local wavelength are used. In ifSTI and related papers, the trigonometric finite 
wave elements (TEWE) is described: the PUM is used with special basis functions 
adapted to waveguides, lasers and geometries with a single dominant wave propa¬ 
gation direction. The finite ray element method of Il79l consists in the use of a PUM 
basis in a first order system of least squares (POSES) formulation; as the unknown 
is constituted by both u and its gradient, more unknowns are needed but the system 
matrix is Hermitian. Pinally, in the hybrid numerical asymptotic method of 1421 . the 
PUM space is constructed by multiplying nodal hnite elements to oscillating func¬ 
tions whose phases are derived from geometrical optics (GO) or geometrical theory 
of diffraction (GTD), e.g. by solving the eikonal equation, cf. il4.2l 

Plane wave bases have been combined in ll97l with the virtual element method 
(VEM) framework mi, in order to design a high-order, conforming method for the 
Helmholtz problem, in the spirit of the PUM, but allowing for general polytopic 
meshes. The main ingredients of the resulting PW-VEM are (i) a low frequency 
space made of low order VEM functions, which do not need to be explicitly com¬ 
puted in the element interiors, (ii) a proper local projection operator onto a high- 
frequency space made of plane waves, and (Hi) an approximate stabilisation term. 
The implementation of the PW-VEM does not require computation of volume in¬ 
tegrals, and no quadrature formulas are required for the assembly of the stiffness 
matrix, for meshes with flat interelement boundaries. 

The hybridizable DG method of employs two discontinuous discrete spaces 
(one scalar and one vector) and a space of Lagrange multipliers on the mesh in¬ 
terfaces. Though Trefftz spaces might be used with this formulation, the authors 
consider basis functions constructed as products of polynomials and geometrical 
optics-based oscillating functions, similar to those in ll42l but discontinuous. 

A Trefftz approach has been proposed in the context of finite difference schemes 
in the flexible local approximation method (ELAME) by I. Tsukerman, see e.g. the 
comprehensive review 111 131 . In the ELAME, the Taylor expansion of the solution 
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to be approximated used to define classical finite difference schemes is substituted 
by an expansion in a series of Trefftz basis functions, leading to better accuracy. 

Oscillatory basis functions have been successfully used in boundary element 
methods, in particular for scattering problems, see the review on the hybrid numer¬ 
ical-asymptotic BEM (HNA-BEM) II20L the plane-wave basis boundary elements 
1961 §3] and the extended isogeometric boundary element method (XIBEM) 0931 . 


3 Trefftz discrete spaces and approximation 

Given a Trefftz variational formulation of a BVP, as those in the definition of a 
Trefftz finite element method is completed by the choice of a discrete space 

Vp{.%) = {v€ T{^,) : V|^ G Vp^{K)] C 

where (X) is a p/f-dimensional space of functions v on X such that 4 v + k^v = 0. 
We describe next the main features of the most common local Trefftz spaces Vpf,{K)-, 
we do not consider Lagrange multiplier spaces on mesh faces for the methods in 
112.2.31 The discussion of the conditioning properties of the basis functions described 
and of the techniques for their numerical integration is postponed to 0 


3.1 Generalised harmonic polynomials (GHPs) 

Generalised harmonic polynomials are smooth Helmholtz solutions that are separa¬ 
ble in polar and spherical coordinates in 2D and 3D, respectively, i.e. circular and 
spherical waves (also called Eourier-Bessel functions or Eourier basis). The local 
spaces Vpf,{K) are defined as follows: 

2D: Vp^{K) = Iv: v{x)= £ aiJi{k\x-XK\)e’^^, Ue €C\, 

e=-‘iK 

3D: yp^(X) = |v: v(x) = ^ ^ |x-x;f|)T/’G c|, 

i=0m=—i 

where xk € K (e.g. is the mass centre of K), 9 is the angle of x in the local polar 
coordinate system centred at xk, Ji is the Bessel function of the first kind and order 
£, is a basis of spherical harmonics of order f (see e.g. If85l eq. (B.30)]), 

and ji is the spherical Bessel function defined by ji(z) = The space 

dimension pK is given by pK — 2qK + 1 in 2D and by pK = {qk + 1 in 3D. We call 
qK, the maximal index of the (spherical) Bessel functions used, the “degree” of the 
GHP space, as it plays the same role of the polynomial degree in the approximation 
theory. A particular feature of GHP spaces is that they are hierarchical. 
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The name “generalised harmonic polynomials” was coined in IISOll and comes 
from the fact that they are images of harmonic polynomials under the operator that 
maps harmonic functions into Helmholtz solutions, in the framework of Vekua- 
Bergman’s theory 0121111411 (see also 05011870 1. The same theory allows to transfer 
approximation results for harmonic functions by spaces of harmonic polynomials 
into results on the approximation of Helmholtz solutions by GHPs. The density of 
GHPs in a space of Helmholtz solutions was proved in ll50l Th. 4.8] and 111 141 §22.8]. 
Approximation estimates in two dimensions were first proved in ll28] Th. 6.2] (in 
norm) and in 1801 (in Sobolev norms), and later sharpened and extended to three 
dimensions in li86l . We summarise here the estimates of 18^ Th. 3.2]. 

Let D G R", n = 2,3, be a bounded, open set with Lipschitz boundary and diame¬ 
ter hf), containing Bpij^{xj)) (the ball centred at some Xjj G D and with radius pho), 
and star-shaped with respect to where 0 < po < p < 1/2. Assume that 

u G (D), s S N, satisfies Au + k^u = 0 in D and define the k-weighted Sobolev 
norm := J ^ N, where is the Sobolev semi¬ 

norm of order m on D. 

i) If n = 2 and D satisfies the exterior cone condition with angle li8^ Def. 3.1] 
(Xo = 1 if D is convex), then for every L> s there exists a GHP Ql of degree at 
most L such that, for every / < s + 1, it holds 

where the constant C > 0 depends only on the shape of D, j and s, but is inde¬ 
pendent of ho, k, L and u. 

ii) If n = 3, there exists a constant Xd > 0 depending only on the shape of D, such 
that for every L > max{s,2*/'^®} there exists a GHP Qt of degree at most L such 
that, for every j < s -f 1, it holds 

where the constant C > 0 depends only on the shape of D, j and s, but is inde¬ 
pendent of hj), k, L and u. 

The main difference between the two results is that the positive shape-dependent 
parameter Xd entering the exponent of L (thus the p-convergence order) is explicitly 
known in 2D (it depends on the largest non-convex corner of D) but not in 3D. 

Exponential convergence of the GHP approximation of Helmholtz solutions that 
possess analytic extension outside D were proved in ll^ Prop. 3.3.3] and improved 
in 2D in li54ll . based upon the corresponding result for harmonic functions of li55ll . 
Roughly speaking, the error is bounded by a negative exponential of the form 
Cexp(—feL) ^ Cexp(—*^), while classical bounds for polynomials achieve 

at most Cexp(—since the dimension po of the GHP space of order L is 
while the dimension po of the polynomial space of degree L is ^(L"). 
Thus, Trefftz methods based on GHPs (and similarly on PWs) can achieve better 
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asymptotic order than standard schemes; however the value of the positive coeffi¬ 
cients b, C and their dependence on the B VP and discretisation are not entirely clear. 

Approximation estimates in the (discontinuous) spaces Vp{^h) immediately fol¬ 
low from the local approximation estimates with D = K, for all K € In case 
of (//'-conforming) partition of unity spaces enriched with GHPs, global estimates 
follow from combining the local estimates with lISTl Th. 2.1]. 

GHPs have been proposed in numerous Trefftz formulations; LS ||8^I107L 
UWVF 123, VTCR EH, hybrid-Trefftz ||99l eq. (62)], direct and indirect single¬ 
element schemes ll22lfT2n . HELS Uni, MPS ifTSira . 


3.2 Plane waves (PWs) 

Plane waves probably constitute the most common choice of Trefftz basis functions. 
In this case, the local space Vpj^{K) is defined by 



(14) 


where j C K", |d^| = 1, are distinct propagation directions. To obtain iso¬ 

tropic approximations, in 2D, uniformly-spaced directions on the unit circle can be 
chosen (i.e. d^ = {cos(27t£/ pk), sm(27t£/ pk)))', in 3D, 01031 and ||94l provide di¬ 
rections that are “almost equally spaced” (see lU §3.4] for a simpler version). In 
these cases, the PW spaces are not hierarchical. However, one of the potential bene¬ 
fits of PW approximations is the possibility to depart from the isotropic case and to 
adapt the basis propagation directions to the specific BVP at hand and to different 
elements, either a priori or a posteriori, see il4.2l 

The linear independence of arbitrary sets of plane waves (and of their traces) 
is proved in Hum]. PW bases whose linear independence does not degenerate for 
small values of kh]^ were introduced in P6l §3.1] in 2D and in If86l §4.1] in 3D (see 
also ||85] §3.4.1]) for analysis purposes. These stable PW bases converge to GHP 
bases in the low-frequency limit 18^ p. 815]. The existence of these stable bases, 
which is instrumental to the derivation of approximation estimates for Helmholtz 
solutions in PW spaces in ll8^ . is guaranteed, provided that the set of directions 
j constitutes a fundamental system for certain harmonic polynomials. In 2D, 
any choice of pK = 2qK + 1 distinct directions, qK being the maximal degree of 
the considered harmonic polynomials, guarantees this property. In 3D, sufficient 
conditions on pK = [qK + 1)^ directions are stated in ll8^ Lemma 4.2]. 

Approximation estimates in PW spaces can be derived from similar bounds for 
GHPs such as those in il3.1l In lISOl Ch. 8], GHPs were approximated by PWs 
by approximating their smooth Herglotz kernel with delta functions, leading to p- 
estimates in 2D, while in |f86l the Jacobi-Anger expansion was used to link PWs 
and GHPs in 2D and 3D. Theorems 5.2 and 5.3 of If86l (see also ll85l §3.5]) show 
that Helmholtz solutions of given Sobolev regularity can be approximated in PW 
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spaces with hp-estimates similar to those shown in l|3.1l for GHPs. For PWs, these 
estimates hold with L = so that plays the role of a “degree” for the consid¬ 
ered PW space. As mentioned, for these bounds to hold in 3D, the PW directions 
have to satisfy some further conditions. A different derivation of /r-approximation 
estimates based on a Taylor argument can be found in lfT9l Th. 3.7]. In ll95l . the 
PW approximation of Helmholtz solutions on the unit disc is analysed in detail, 
together with the conditioning of different linear systems used for its computation 
(least squares and collocation for a Dirichlet problem on the disc) and the implica¬ 
tions on the accuracy of the approximation computed in finite-precision arithmetic. 
We refer again to ll54l §5.2] for the exponential convergence in 2D of PW approx¬ 
imations of analytic Helmholtz solutions (see also ll85l Rem. 3.5.8] which holds in 
2D and 3D). 

Similar to PWs are the evanescent waves: the basis elements have the same ex¬ 
pression v(x) = e'^'*** but with a more general dGC",d-d = l.Ifd = d/;-|- id/, with 
d//,d/ S K", then v oscillates in the direction du (with wavenumber k\dR\ > k) and 
decays exponentially in the orthogonal direction d/ (i.e. |v(x)| = *). Evanes¬ 

cent waves are used in combination with plane waves to approximate interface prob¬ 
lems in the DEM El] and the UWVE |77l, and to represent outgoing waves in a 
2D unbounded half-strip of the form {a <x < b,y > c} in E1[II1. 

A special combination of propagative and evanescent waves is typically used in 
the WBM. We describe a 2D version of this space as in ll24l eq. (14)-(21)] (see 
lIZTl §4.1] for 3D). This space is not invariant under rotation but depends on the 
choice of the Cartesian axes. Eor a mesh element K, we fix a truncation parameter 
N >0 (typically 1 < < 6) and define := ^^V{xuy\),(s 2 ,yi)cK l-^i -^ 2 \ and Ly := 

sup^^j (x 2 y 2 )eK bi ~y2\ as the edge lengths of the smallest rectangle containing 
K and aligned to the Cartesian axes. Two sets of basis functions are used: 

cos(;t,,■x)e*‘^/^^ 

e±l^<: -kyjX 

for a total dimension pK = A- + 2{\NkLx/n\ + \NkLy/n\). Each basis function is 
half the sum of two plane (or evanescent) waves, symmetric to one another with 

respect to the x or y axis: e.g. cos{kxjx) exp(i \Jk^ — ^xjy) ~ 

with d^- := {±kxj/k,^l — (kxj/kY). A maximum of A + 2{\kLx/K\ -\- \kLy/K\) 
basis functions are propagative PWs, this number designed to keep the conditioning 
under control. If A > 1, then roughly a fraction {N — \ )/N of the total basis func¬ 
tions are evanescent waves decaying in a direction parallel to one of the Cartesian 
axes. Refinement is obtained by increasing N: for A < 1 only propagative waves are 
present, for higher values evanescent waves are introduced. 

In 2D, both evanescent and plane waves may be written as exp{|(i(v -I-1 /v)x-\- 
(v — 1 /v)y} = exp{iA:(xsin 0 -fycos 0}, parametrised by v G C or 0 g C with v = 


kxj:=jj, 7 = 0,..., [AkLf/:/rJ, 
kyj--=^, 7 = 0 ,..., [NkL’^/n \, 









28 


Ralf Hiptmair, Andrea Moiola and Ilaria Perugia 


e‘®; these waves constitute the test space (but usually not the trial) for the Fokas 
method in II23II105I and ll^ §7.3.4] (see also 112.41 1. 


3.3 Fundamental solutions and multipoles 

Fundamental solutions and multipoles are Helmholtz solution in the complement 
of a point and satisfy Sommerfeld radiation condition (lim^-i-ocrT — iku) = 0, 
where r = |x|). They are particularly useful to dehne Trefftz spaces on unbounded 
elements, e.g. for scattering problems. 

If the local spaces are spanned by fundamental solutions, simple sources are 
located at distinct poles in the complement of K: 

2D : Vpi^{K) = lv: v(x) = ^ a<<H^^'(k|x-Xf|), s c|, 

^ i=\ ’ 

3D: Vp^{K) = \v: —p, 


where is the Hankel function of the first kind and of order 0. Different a priori 
or a posteriori strategies are used to fix the location of the poles, see 112.1. H and the 
references cited therein. As the distance of the points x^ from K increases, these ba¬ 
sis functions approach plane waves, so they permit flexibility not only in the choice 
of the propagation directions but also in the wavefront curvature. 

Apart from the MFS and its modifications (see 112.1. H and lfTll9] [T0][3ni921ll20l ). 
spaces of fundamental solutions have been used in connection to the UWVF (see 
158), where ray-tracing is used to determine the poles, and ll57l '). 

Theorem 6 of 01041 ensures that Helmholtz solutions in K can be approximated 
in Holder norms by fundamental solutions centred at any “embracing boundary” in 
2D and 3D, under weak assumptions on the regularity of dK. We are not aware of 
any result providing orders of convergence. 

An alternative approach consists in choosing local spaces generated by multipole 
expansions, where multiple sources with increasing order are located at a single pole 
xo (or only at few poles): 

2D:Vp^{K) = \v: v{^)= £ a^//;'^(k|x-xo|)e'^®, G c|, 

t=-qK 

3D: Vpf.{K) = \v: v(x) = £ £ {k\s.-XQ\)Y^(y —afo,GC|, 

^|x-Xo|7 J 

where (PP) are Hankel functions (spherical Hankel functions, respectively) of 
the first kind and order 1. As for the GHPs in 113.11 0 is the angle of x in the local 
coordinate system centred at xq, which is located in the complement of K, and the 
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space dimension is = 2qK + 1 in 2D and = {qK + 1)^ in 3D. According to IfTOl 
Rem. 2.2], fundamental solutions lead to more stable methods than multipoles. 

Multipole spaces have been used in connection to LS schemes 09011 1071 . WBM 
El eq. (23)1, EH §4.1.2], hybrid-Trefftz JH eq. (63)], HELS Oil, source simula¬ 
tion techniques 19^ . null-field fTSll and single-element schemes ll8] l^fT2Tll . In 1491 
and related papers, some 2D multipoles with suitably chosen index i (not necessar¬ 
ily integer) are used on infinite sectors, in such a way to ensure continuity of discrete 
functions across rays; this might be more efficient than full multipole spaces for so¬ 
lutions with a preferred propagation direction. 


3.4 Other basis functions 

Other discrete Trefftz spaces have been proposed in literature for use with the vari¬ 
ous approaches covered in 

In 2D, corner waves such as 7^/a(A;|x|)sin(f0/a), with f S N and 0 < a < 2, 
capture the behaviour of Helmholtz solutions near a domain corner of angle na. 
They have been used e.g. in the WBM l25l . in LS methods UlOlllOTiri 191 and in the 
MBS |[T6l[36l. In 11201 . they are used with a = 2 on tips of ID screens to repre¬ 
sent the strong singularities of the solution in a non-Lipschitz domain. Theorem 6.3 
of l28l uses Vekua-Bergman theory to give orders of convergence for the approxi¬ 
mation of singular functions by spaces of corner waves and GHPs (see also ifTol §5] 
and references therein). We are not aware of any use of similar functions in 3D. 

The wave band functions, introduced in the VTCR context 11011 . are Herglotz 
functions with piecewise-constant kernel, e.g. 2D. 

In the presence of a circular hole, suitable combinations of Hankel and Bessel 
functions a priori fulfil homogeneous boundary conditions 11071 eq. (13)]. 

If the wavenumber varies inside an element, the basis functions described so 
far do not lead to Trefftz methods. In case of linearly variable wavenumber. Airy 
functions can be used to construct Trefftz spaces 11 101 . In 1641651 generalised plane 
waves in the form for suitable polynomials P, are introduced and analysed in 
a UWVF setting; they solve a perturbed Helmholtz problem and converge with high 
orders in Iik- Similar “almost-Trefftz” waves are used in l43l and named oscillated 
polynomials. Modulated plane waves, i.e. products of PWs and polynomials, are 
the basis functions of the DG method of I141I15I : as they are only “approximately 
Trefftz”, volume terms appear in the formulation. 

Products of (continuous) low-order polynomials and PWs or GHPs constitute the 
basis of the PUM 151I73|[8TI|9^108I . while products of polynomials and oscillating 
functions derived from high-frequency asymptotics are basis elements in 142119 1 1 . 
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4 Further topics 

4.1 Assembly of linear systems 

All the Trefftz finite element methods for O discussed in ij2] give rise to dense 
or sparse linear systems of equations. Entries of coefficient matrices are obtained 
by integrating products of (derivatives of) trial and test functions over bounded d- 
dimensional sub-manifolds of Q,d <n. The stable and accurate (approximate) eval¬ 
uation of these integrals is a key implementation issue. 

Among all Trefftz approximation spaces and associated bases presented in ^ 
plane waves (PWs) * (either propagative with d € M" or evanescent with d S 
C") are exceptional, because they allow a closed-form evaluation of their integrals 
over any flat sub-manifold with piecewise flat/straight boundary. For instance, in all 
variants of PW-based Trefftz methods on polyhedral meshes in 3D, expressing mesh 
faces by 2D parametrisations, we eventually encounter integrals of the form 



E C a bounded polygon, w e constant. (15) 


Then we can take the cue from ll3^ §2.1] or ||29] §4] and apply integration by parts 
in order to reduce (flSl l to integrals over the straight edges e\,e 2 , ■ ■ - Cq, q & of F: 



where is the exterior normal at e^. Then, as in ii ch. 2], if ei = [a,b], 
a,b e R^, we find, /g^exp(w-x)ds = exp(w-a)|b —a|v/(w- (b —a)), where v/(z) := 
(exp(z) — l)/z. Of course, a numerically stable implementation of this function for 
small arguments is essentiaO This approach can be generalised to yield analytic for¬ 
mulas for computing integrals of products of PWs times polynomials, see II29II38I . 
with increased computational effort, however. 

Approximate evaluation of the integrals becomes inevitable for all choices of Tr¬ 
efftz basis functions other than PWs, and even for a PW basis on meshes with curved 
elements. Then Gauss-Legendre numerical quadrature seems to be the most widely 
used option. However, the integrands may be oscillatory, which delays the onset of 
(exponential) convergence of the quadrature error until the number of quadrature 
points surpasses a threshold roughly proportional to the ratio of the local mesh size 
and the wavelength. This leads to higher computational cost per degree of freedom 
for larger values of khK- One may think of using special quadrature rules for os¬ 
cillatory integrals, as derived, for instance, in ll62l . Those avoid an increase in the 
number of quadrature points for growing spatial frequency of the oscillations, but 
unfortunately require precise knowledge of the oscillatory term in the integrand. 


* A stable algorithm for point evaluations of y/ even for arguments close to 0 is provided by the 
MATLAB function expml. 
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Besides classical h-, p- or /ip-adaptivity, Trefftz methods offer scope for more so¬ 
phisticated adaptive strategies consisting in the choice of specihc basis functions for 
different BVPs and in different mesh elements, either a priori or a posteriori. 

The main strand of a priori adaptive Trefftz methods falls into the category of 
hybrid numerical-asymptotic methods. High-frequency limit models, such as ray 
optics or geometric theory of diffraction (GTD), guide the selection of local Trefftz 
spaces in the individual cells of a mesh. In a non-Trefftz PUM framework this idea 
was pursued in Il42ll . and within the hybridizable DG method in ||9T], in both cases 
for 2D acoustic scattering at a smooth sound-soft object. In these works, local phase 
factors X exp(ik5(x)) derived from reflected and diffracted waves multiply stan¬ 
dard continuous nodal basis functions, in ll42l . or local polynomials, in 11911 , thus 
generating a basis for (local) trial spaces. 

The policy of incorporating local directions of rays is particularly attractive for 
PW-based methods, because PW basis functions naturally encode a direction of 
propagation. For problems where excitation is due to an incident PW and mate¬ 
rial properties are piecewise constant, ray tracing and related techniques 1911 §3.2] 
based on geometric optics (specular reflection and Snell’s law of refraction at ma¬ 
terial interfaces) can provide information about the local orientation of wave fronts 
for k-^-oo. PWs matching the found ray directions are then used to build local bases, 
either exclusively or augmented by a reduced set of generalised harmonic polyno¬ 
mials (GHPs) or “equi-spaced” PWs. 

This idea for TDG was first outlined and tested in Ql and further elaborated 
and extended in ll58] Ch. 5] (for UWVF). In the latter work, in an attempt to resolve 
curved wave fronts and take into account diffracted waves from corners, also Han- 
kel functions x i—!> //Q*^(k|x — y*|) with y* outside a mesh cell have been proposed 
as local basis functions. Approximation of curved wave fronts deduced from GTD 
corrections is also attempted in El. There the authors move beyond Trefftz meth¬ 
ods and use DG with trial spaces of polynomially modulated PWs, which are more 
suitable for approximating propagating circular waves. 

In simple 2D situations with convex smooth or polygonal scatterers and incident 
plane wave, overall accuracy seems to benefit substantially from a priori directional 
adaptivity. However, if there are more than only a few dominant wave directions 
as in the case of more complicated geometries, trapping of waves, dark zones and 
shadow boundaries, current directional adaptivity soon meets its limitations. On the 
other hand, this strategy appears as the most promising way to achieve k-uniform 
accuracy with numbers of degrees of freedom that remain k-uniformly bounded 
or display only moderate growth as k —>• oo. The potential of this idea has been 
strikingly demonstrated in the case of BEM for 2D scattering problems ll20l . 

A posteriori directional adaptivity seeks to extract information about dominant 
wave directions from intermediate approximations of u. A reflne-and-coarsen strat¬ 
egy is embraced in IflTll . In each step of the adaptive cycle it first computes a PWDG 
solution u of the scattering problem based on a relatively large number of local Tr¬ 
efftz basis functions (GHPs and PWs). Subsequently, by solving local non-linear 
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L^-least squares problems, the directions of fewer PWs are determined so that u can 
still be well approximated locally. 

A p-hierarchical error indicator is studied in ll44l . In a step of the adaptive scheme 
starting from the approximate solution u a presumably improved solution u is com¬ 
puted using double the number of local PWs. Then a single local plane wave direc¬ 
tion ds: on a mesh element K is extracted from the error e(x) := m(x) — m(x) through 
the projection formula 



Detailed numerical experiments are reported in ll44l Ch. 6]. In the pre-asymptot- 
ic regime, when the resolution of the trial spaces is still rather low, one observes 
a pronounced gain in accuracy in the case of the adaptive approach compared to 
approximation with the same total number of equi-spaced PWs. 

Directional adaptivity for Trefftz methods has also been tried in other flavours. 
In the context of least squares methods as discussed in 112. II an offset angle for the 
sets of local equi-spaced PWs is introduced as another degree of freedom in 13 , 
aiming to align them with a local dominant wave direction. For the VTCR method 
presented in 112.3.11 an error indicator based on local wave energy is used in 11021 
to steer angular refinement of local Trefftz spaces. 

A posteriori mesh adaptivity is considered in ll66l . where classical “elliptic” 
error estimation and mesh rehnement strategies are adapted for the /z-version of 
TDG. In a low-frequency setting, the method inherits the good performance of the 
underlying adaptive mesh refinement algorithms for polynomial DG for the Poisson 
equation. However, there is little hope that this carries over to larger wavenumbers k. 
A similar error estimator, aimed at adaptive mesh refinement, has been described 
in 0 §3.2] for the DEM/DGM presented in ^Z23] 


4.3 Ill-conditioning and solvers 

The linear systems of equations spawned by PW-based finite element methods are 
highly prone to ill-conditioning, when high resolution trial spaces are used, see 
e.g. ED §5], ED §4 .3], iBOl, and ||72l for a PUM setting. This is largely caused by 
an inherent instability of the PW basis on cells, whose size is relatively small com¬ 
pared to the wavelength. Intuitively, for |x| <C the functions x i—e*^**'"* from 
d) are almost constant, hence, nearly linearly dependent, cf. 11721 §4.2]. The same 
heuristics applies, when their density increases; even for cell sizes comparable to the 
wavelength, PWs are hardly distinct when their directions are close, cf. m §4.3]. 

Empirically, for the local PW Galerkin matrix M/f associated with the if inner 
product on a single mesh cell K, we find that its spectral condition number grows 
like ^ for cell size Hk 0, where q > 0 is proportional to the number pK 
of (approximately uniformly spaced) PWs in 2D, and to the square root of pK in 
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3D. Essentially, q is related to the “degree” of the considered set of PWs; see 
113.21 Even worse, the condition number soars exponentially in q: cond 2 (M; 5 -) ~ e“^ 
for q oo and a > 0; see Appendix A. A similar explosion of condition numbers 
is observed for the full systems matrices as meshes are refined or more PW basis 
functions per element are used. 

There is circumstantial evidence that direct sparse elimination can cope fairly 
well with the ill-conditioned linear systems arising from UWVE or PUM, see Bol 
§5.3.3], Gil. Yet, eventually the instability of the basis will impact the quality of 
the solution 111081 §5.4]. A remedy proposed in for the UWVE is to limit pK 
based on monitoring condition numbers of element matrices. Apparently, this also 
curbs the condition number of the global system matrix. Alternatively, there exist 
different heuristic recipes for choosing a priori the number of PWs per element to 
balance accuracy and conditioning; in 2D, the widely cited ll60] eq. (14)] suggests 
PK = m\md{khK + C{khKY/^) with 3 < C < 14 for the UWVE, while ET] §5.1.1] 
proposes pK = \ 2khK\ for the VTCR. Eor the WBM, li24l §3.2] proposes a rule to 
balance propagative and evanescent basis functions, see il3.2l 

The most straightforward cure for instability would trade the PW basis of {K) 
from (fl4l i for a more stable basis, found by local orthonormalisation as in the case of 
polynomial EEM, cf. the approach from li^ §3.1]. However, instability may sneak 
in through the back door and manifest itself in severe impact of round-off errors 
during orthonormalisation and recombination of element matrices. The use of high- 
precision arithmetic may be advisable, but has never been documented. 

Eor the sake of stability, PWs may be replaced by the generalised harmonics 
polynomials introduced in 113.11 In 2D, a scaling of the GHPs has been devised 
in llTTl , in order to lower the condition number of the resulting UWVE: 

k^\j[{khK)\^ + \Jt{khKf 

In Cll, it is also shown that the conditioning of GHP-based UWVE schemes is 
better than for methods based on PWs, and that it improves on regular meshes. This 
might be related to the orthogonality of GHPs on balls. 

The numerical experiments in li58] §3.7] suggest that the use of fundamental 
solutions as basis functions may considerably reduce the conditioning of UWVE 
matrices, at the expense of accuracy. Both accuracy and conditioning increase the 
further the centres of the fundamental solutions are from the element. 

The use of iterative solvers for linear systems generated by Trefftz methods en¬ 
tails preconditioning. Eor PW basis functions, the first proposal in lfT9] §2.4] for the 
UWVE was a local preconditioner, equivalent to an orthonormalisation of the PW 
basis with respect to an inner product on the boundary of mesh cells. An interest¬ 
ing connection of the local preconditioner with non-overlapping optimised Schwarz 
domain decomposition methods was discovered in na. The local preconditioner 
was used in conjunction with a BiCGStab Krylov subspace solver in M and aug¬ 
mented by a coarse-grid correction in the spirit of non-overlapping domain decom¬ 
position in M59II1181 . The coarse space is again spanned by PWs. This is also true for 
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the two-level sub-structuring preconditioner proposed for DEM/DGM (see il2.2.3l l 
in 1341 . Two-level, non-overlapping Schwarz domain decomposition precondition¬ 
ers for PWDG (essentially UWVF) have been tested in Q; these preconditioners 
seem to be robust with respect to the wavenumber k and the local number of PW 
directions, although they do not seem to be perfectly scalable with respect to the 
number of subdomains. 


5 Assessment and conclusion 

Faced with a flurry of different Trefftz methods and a wealth of numerical data, we 
feel at a loss about making unequivocal statements about the merits of Trefftz meth¬ 
ods, let alone ranking them according to some undisputed criteria. Rigorous theory 
is available for LS methods ( 112. II) . TDG ( 112.2.11) . and PUM ( 112.51 ). Combined with 
approximation results for suitable Trefftz bases, this leads to better asymptotic esti¬ 
mates in terms of orders of convergence in the number of degrees of freedom to what 
is available for polynomial FFM (e.g. II52II54II ). The dependence of crucial constants 
on the wavenumber k is explicit in several cases, but the orders in k are usually not 
better than for polynomial methods. Thus theory fails to provide information about 
the key issue of “k-robust” accuracy with “k-independent” cost. Moreover, numeri¬ 
cal dispersion will also haunt local Trefftz methods in the case of /r-reflnement; thus 
they provide no escape from the pollution error. 

We also advise caution when reading numerical experiments, because they may 
be tarnished by selection bias, making authors subliminally pick test cases matching 
the intended message of an article. Disregarding this, even “objective” comparisons 
are inevitably confined to a few simple model problems. This is problematic, be¬ 
cause different model problems sometimes seem to support opposite conclusions. 

From our experience, the power of Trefftz methods can best harnessed by p- 
reflnement using approximation by Trefftz functions in regions as large as possi¬ 
ble. In the presence of singularities we recommend either the use of corner basis 
functions ( 113.41 ) in 2D, or kp-reflnement, maybe using standard polynomial approx¬ 
imation on small elements as in ISSlI . There is a solid theoretical foundation, when 
this is done in the LS, TDG, or PUM framework. The resulting methods should be 
able to compete successfully with polynomial FFM even in their more sophisticated 
versions tailored to wave propagation problems 0301135118^ . 

The discussion of adaptive approaches in li4.2l hints that some Trefftz trial spaces 
have approximation capabilities well beyond the reach of polynomials. Directional 
adaptivity seems to be very promising, but much research will still be required to 
convert it into a reliable practical algorithm. The same applies to iterative solvers 
and preconditioners for Trefftz schemes, see 114.31 which might also benefit con¬ 
siderably from the extra information contained in Trefftz trial spaces. Hence, we 
believe that many exciting possibilities offered by the idea of Trefftz approximation 
still await discovery and that the full potential of Trefftz methods is only gradually 
being realised. 
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Appendix A: Condition numbers of plane wave mass matrices 

Given a wave number k>Q and p gN distinct unit vectors C K", £ = l,...,p, and 
a domain K cM" with barycentre x^, the symmetric positive definite plane wave 
element mass matrix on K is defined as 



For n = 2 we computed spectral condition numbers of Mjf for equi-spaced directions 
= (cos(2;r£/ p),sm{2n£/p)), i = 0,... ,p — 1. For n = 3 we choose the directions 
d^ as the “minimum norm points” according to I.H. Sloan and R.S. Womersley 
II103II116l . These points are indexed by a level ^ G N and p = (g'+ 1)^. The spectral 
condition numbers are plotted in Figure [T] for n = 2, K = (—1,1)^, and Figure |2] 
for n = 3, K = (—1,1)^. They have been computed with MATLAB using the high- 
precision arithmetic (200 decimal digits) provided by the Advanpix Multi-Precision 


Toolbo?(0. 


^ http://www.advanpix.com/ 
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Fig. 1: Condition numbers of element mass matrices on the square (—1,1)^ 




Fig. 2: Condition numbers of element mass matrices on the cube (—1,1)^. 








































